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Abstract 

Atomistic simulations are used to study the bending of rectangular graphene nano ribbons 
subjected to axial stress both for free boundary and supported boundary conditions. The shape 
of the deformations of the buckled graphene nano ribbons, for small values of the stress, are sine 
waves where the number of nodal lines depend on the longitudinal size of the system and the 
applied boundary condition. The buckling strain for the supported boundary condition is found 
to be independent of the longitudinal size and estimated to be 0.86%. From a calculation of the 
free energy at finite temperature we find that the equilibrium projected two-dimensional area of 
the graphene nano ribbon is less than the area of a flat sheet. At the optimum length the boundary 
strain for the supported boundary condition is 0.48%. 

1 Introduction 

Graphene, is a newly discovered almost flat one-atom-thick layer of carbon atoms which exhibits 
unique electronic properties and unusual mechanical properties [1, 2]. Recent experiments showed 
that compressed rectangular monolayer graphene on a substrate with size 30 x 100 fim 2 is buckled 
at about 0.7% strain [3]. Moreover tensional strain in monolayer graphene affects the electronic 
properties of graphene. The strain can generate a bulk spectral gap in the absence of electron-electron 
interactions as was found within linear elasticity theory, and a tight-binding approach [4]. Different 
morphological patterns of carbon nano-structures subjected to external stress were obtained by using 
atomistic simulations [5]. Furthermore, atomistic simulations showed that the Young modulus and 
the fracture strength decrease only weakly with the width of the graphene nano ribbon [6]. 

In this paper we study the deformations and the stability of rectangular monolayer graphene nano 
ribbons (GNR) subjected to axial stress using atomistic simulations and the Jarzynski theorem to 
calculate the free energy [7] . Recently, Colonna et al applied the free energy integration based method 
to explain the melting properties of graphite [8]. We will compare the obtained critical buckling force 
with the one predicted by elasticity theory. We found several longitudinal deformation modes and 
predict that the axial buckling boundary strain is independent of the size in the case of laterally 
supported GNRs. Moreover, from a calculation of the free energy, uncompressed GNR is thermody- 
namically less stable than the GNR at the buckling threshold. But the buckled state is less stable 
than the GNR at its optimum length. 

This paper is organized as follows. In Sec. 2 we introduce the atomistic model and, the simulation 
method. Section 3 contains a discussion of the elasticity theory predictions for both free boundary 
condition and laterally supported boundary condition. We also give the buckling thresholds and 
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Figure 1: (Color online) Schematic model for a plate under axial strain with free boundary condition 
(a model for elasticity theory), dashed-dotted rectangle is the initial non-compressed plane. We list 
all relevant variables describing the GNR. 



the obtained deformations and compare our results to those from elasticity theory. Results for the 
Young's modulus and pre-stresses are presented and compared to available experimental results. The 
stability of buckled GNRs are studied in the last part of Sec. 3. In Sec. 4 we will conclude the paper. 

2 Method and model 

Classical atomistic molecular dynamics simulations (MD) is employed to simulate compressed GNRs 
using Brenner's bond-order potential [9]. Our system is a rectangular GNR with dimensions a x 6, 
in x and y directions with armchair and zig-zag edges, respectively. For simplicity we set the lateral 
dimension b = 10n y y/3ao where ao=0.142 nm, n y = 3,4 and the longitudinal dimension a = 30n x ao 
with n x — 2,3, ...10. In Fig. 1 we depict a schematic model for a GNR under axial strain with free 
boundary condition (at y = and y = b) and list all relevant variables describing GNR under axial 
strain. Note that n x and n y are two integer numbers that are related to the number of atoms in the 
armchair and zig-zag direction, respectively. The corresponding values for length and width of GNRs 
can be found in Table 1. 

Initially the coordinates of all atoms are put in a flat surface of a honey comb lattice with nearest neighbor 
distance equal to ao and the initial velocities were extracted from a Maxwell-Boltzman distribution at the given 
temperature. Before starting the compression, the system is equilibrated during 50 ps (100.000 time steps). 
Compressing direction is always x and two rows of atoms in both right and left edges, x = and x — (2, are 
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Table 1: Length a and width b of GNRs which are related to the integer numbers n x and n y through 
b — 10n y y/3ao and a — 30n x ao where ao=0.142 nm. 
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Figure 2: (Color online) (a) Buckled graphene nano-ribbons for four different length n x . Here 
Ax =0.768 nm. (b) Variation of the out of plane displacement in log-scale for GNRs with two 
different lengths versus strain, (c) Symbols are the calculated buckling force in log-log scale and 
dashed line is the curve f c 



7T 2 P 

a 2 



fixed during the compression steps (Sx = 0.04 A) with the rate fi = 1.6 m/s. The boundary axial strain after 
t compression steps is 

e x = te , e = 2Sx/ (30a n x ), (1) 
where eo = is the strain after a single compression step. After each compression step, we wait 2.5 ps 

Tlx 

to allow the system to relax. For the edges at y = and y = 6, we used the supported boundary condition 
and the free boundary condition. We simulated the system at room temperature and employed a Nos'e-Hoover 
thermostat. 



3 Buckling graphene nano-ribbons and comparison with elasticity 
theory 

3.1 Free boundary condition 

For a simple bar with length a, under an axial symmetric load applied at its ends, classical Euler's column 
equation describes basically the buckling problem [10]. Governing differential equation for the deflection value, 
5w, becomes the harmonic oscillator equation 

6w" + n 2 5w = 0. (2) 

Here n 2 = f c /P where f c is the buckling force (or critical force), P is a parameter which is related to the Young's 
modulus and the moment of inertia of the rod cross-sectional axis that is perpendicular to the buckling plane. 
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The solution is Sw = Asm(tzx) + Bcos(tzx). For the boundary condition with zero deflection at the ends, we 
have 5w = Asin( 2 ^ L x) which are sine waves. Substituting k = ^ yields the buckling force, 

n 2 ir 2 P 

fc = ^^ (n = l,2,...). (3) 

The buckling stress can be written as <j c = |, where S is the area of the cross section of the bar. If the bar is 
thin enough or long enough buckling can happen elastically independent of the type of material [10]. 

The shape of the lowest mode (n=l), after and close to the buckling threshold, is a half sine wave. Higher 
modes are possible only if the column is physically constrained from buckling into the lower modes by supporting 
mid points [10]. 

For GNRs under axial compression with free boundary condition we focused on the system with n y =3. 
After many compression steps GNRs starts to buckle, but the shape of the deformed GNRs depends on their 
size. Figure 2(a) shows some snap shots for the deformed GNRs with various sizes when Ax = 0.768 nm, 
i.e. beyond the buckling threshold. We will discuss on the obtained shapes later. By measuring the buckling 
threshold, we find the forces which cause a sudden change in the shape of GNRs. This threshold is found 
from a direct visualization of the nano-ribbon, and in addition from the sudden increase in the average out of 
plane displacement of the GNR atoms ((h 2 )). The variation of (h 2 ) versus e x for systems with n x = 3, and 
10 are shown in Fig. 2(b). The longer GNRs has a smaller buckling strain, e c ~ 1.2%, 2.1% for n x = 10, 3, 
respectively. 

On the other hand, as can be seen from the simulation snap shots in Fig. 2(a), even for large samples, 
in contrast to the large deformations along the x direction, the deformations along the y direction for each x 
value and also deformations at the boundaries (y=0 and y = b) for each x value are negligible. Therefore the 
rod assumption for GNRs under axial strain is a good approximation. Now, considering the GNR as a rod 
with length a and estimating the buckling forces, allow us to calculate the variation f c versus GNR length, i.e, 
a. Furthermore, note that throughout the present paper, we calculate the force per width. Since for n x < 5 we 
found only deformations with n=l in the beginning of the buckling, hence using Eq. (3) is justified. Fig. 2(c) 
shows the variation of the critical load versus a in log-log scale where the dots are from our simulation results 
and the dashed line is 1L ^~- Here we used Eq. (3) for the buckling force, but, since our system is not a rod, the 
definition of an axial moment of inertia is not meaningful. Therefore, we define an effective bending constant 
for GNRs, i.e. P. Taking the parameter P as a fitting parameter, we found P ~ 0.56 eV. The corresponding 
number for a covalent carbon bond is 4.0 eV [11]. Clearly the bending stiffness for a single covalent C-C bond 
should be much larger than a long rod made of GNR. 

Here we discuss the shapes of the deformation. As we mentioned above, at the start of the buckling the 
shape of systems with n x < 5, is found to be almost the same as for a buckled rod weith n=l, which is similar 
to half sine waves. Note that we should neglect the few rows of atoms at both ends which are fixed during the 
compression resulting in flat ends (see section II). By continuing the compression up to the time 0.5 ns, GNRs 
acquire a parabolic shape. The top two pictures in Fig. 2(a) show two snap shots of the obtained deformations 
after the buckling threshold when e x = 2 x = 6%, 3.6% for n x = 3, 5, respectively. For a fixed reduction in 
the length, Ax, the system with larger n x has a smaller amplitude, i.e. 2.609 nm and 3.412 nm for n x =5 and 
n x =3, respectively. For larger GNRs (n x > 6), i.e. a is larger than 21.3 nm, the deformations (two bottom 
figures in Fig. 2(a)) firstly show the next higher mode with n ~ 2 and after a longer simulation time (depending 
on the size) they transit slowly to the mode n = 1. When the length of GNRs exceeds almost 22 nm, half of this 
size is comparable to the chrematistic length of graphene, i.e. 8-10 nm [12]. Since the chrematistic length is a 
measure of the range over which deformations in one region of graphene are correlated with those in another 
region, so the applied boundary stresses on the edges do not affect the system beyond this characteristic length 
and in the beginning of the buckling we expect the n ~ 2 mode. For GNR with n x > 6 we did not find a simple 
relation between the critical buckling force and the length. For larger systems (n x > 10) the deformations are 
no longer sine waves at least during our simulations time. 

Before ending this section, we calculate the stress-strain curve. Before the buckling threshold, we found a 
linear relation between stress and strain, i.e. <j x = Ee x + <Jo, where ctq is the pre-stress in the system [1]. The 
linear relation is valid for small strains [1, 13]. In our simulations, when the GNRs are not flat and thus not 
compressed, they are not in equilibrium and some boundary tension exist, i.e. pre-stress [1]. For instance, when 
n x = 8 we calculate the applied stress on the right hand side (RHS) edge using ^J- (by using the thickness of 
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Figure 3: (Color online) Stress-strain curve for the GNR with n x —S under axial stress before the 
buckling threshold where. 



graphene equal to h$={).l nm [14]) and we show the obtained stress-strain curve in Fig. 3. The dashed line is 
the fitted line. The slope of this line gives us Young's modulus 1.3±0.07 TPa and ao=-3.3±0.2 GPa (negative 
sign indicates the direction of compression, i.e. -x). For the other GNRs we found Young's modulus in the 
same range (e.g. E = 1.1±0.08 TPa and <Jo=-3.3±0.2 GPa for n x = 10 etc). These numbers are comparable 
to those found experimentally [1]. 

3.2 Supported boundary condition 

For a rectangular plate subjected to the supported boundary condition (when movements at x=0, x = a, y =0 
and y = b along y and z directions are not allowed), elasticity theory [10] tells us that the governing equation 
for the buckled rectangular under uniform compressive axial load per width (/) in the ^-direction, can be 
written as 

DV 4 (Sw) + f(Sw xx ) = 0, (4) 

where 8w is the transverse deflection, 5w xx is the corresponding curvature and D = Eh^/ (12(1 — v 2 )) is the 
flexural rigidity of the plate with thickness ho and Young's modulus E. The general solution for the deflection 
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Figure 4: (Color online) (a) Buckled graphene nano-ribbons for four different length subjected to a 
lateral supported boundary condition. Here applied strains are larger than the buckling strains, (b) 
Variation of the out of plane displacement in log-log scale for GNRs for three different lengths versus 
strain, (c) Force applied on the boundaries for three different sizes: n x — 2,4,6. Horizontal arrow 
show the instability point and vertical dashed line shows the buckling point. 



can be written as a double Fourier series 

oo 

5w = Wm n sm(n7Tx/a) sin(ra7n//6), (5) 

m,n=l 

where (m, n) are integers in order to satisfy the supported boundary condition and w mn is the amplitude of 
each mode (ra,n). Including the appropriate strain energy and using Eq. (5), buckling occurs when [10] 

f mn = ^-[^) 2 + (j) 2 } 2 - (6) 

Lowest value of f mn with respect to the two discrete variables m,n gives the buckling force, f c . The minimum 
buckling force for the considered systems always occurs for m = 1 and various values of n. It is equivalent to 
a single half wave in the lateral ^-direction and various harmonics n in the compression direction, i.e., x. 

We performed several atomistic simulations for GNRs under supported boundary condition for different 
n x and fixed n y =3. In Fig. 4(a), we depict four typical snap shots for different sizes of the buckled GNRs 
(where e x is always larger than the buckling strain, e c ). As can be seen from this figure, those satisfy the 
condition m=l, and by increasing the size of the system, we obtained higher longitudinal modes. Furthermore, 
higher strains (e x > e c ) increase the amplitude of the deformations and increase the number n slightly (usually 
n' = 71 + 1,71 + 2). Similar to the free boundary condition case, the buckling thresholds, for smaller n Xl are 
obtained at smaller axial strains. In Table 2, we list our calculated buckling forces per width (/ c ), buckling 
strains and the prediction from elasticity theory according to Eq. (6). To calculate / mn , we used E = 340 
N/m 2 and v = 0.165 [1] with ho ~ 0.1 nm, as a typical thickness for GNR. Larger thickness yields a larger f mn 
so that ho = ao yields a better agreement between f' c and f mn . Note that elasticity theory is not applicable 
to GNRs under strain [3] for the critical values for force or strains. For instance, although elasticity theory 
predicts that the very small thickness of graphene yields a zero flexural rigidity but the bond-angle effects on 
the interatomic interactions of graphene (three body terms in Brenner's potential) gives a non zero flexural 
rigidity for graphene [3] (see Refs. [28,29] in Ref. [3]). Therefore, because of the non-zero flexural rigidity in 
graphene, larger critical forces with respect to the predictions from elasticity theory, are to be expected. 

Figure 4(b) shows the variation of (h 2 ) versus e x for systems with n x = 3,6,10 when n y = 3. Vertical 
dashed line shows the transition points to the buckled state. Clearly the behavior of (h 2 ) is the same for all 
cases. As we see, surprisingly, the buckling strain is independent of the longitudinal size (vertical dashed line 
in Fig. 4(b)) of GNRs and it varies in the range [0.84-0.89]% while for larger width, e.g. n y = 4, the buckling 
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Table 2: The periodicity number (n) of the sine waves observed in the buckled laterally supported 
GNRs in the longitudinal direction versus the length of the GNR. The calculated buckling force from 
our MD calculations (/ c ) and results from elasticity theory according to Eq. (6), / mn , for two widths 
n y — 3, 4. When calculating f mn we used ho = 0.1 nm and ho = ao was used for f' mn . e c is the strain 
at which the GNR buckles. 

strain varies in the range [0.8-0.93]%. The average for n y = 3 is 0.8688% and for n y = 4 is 0.8677% which are 
very close. Therefore we conclude that the buckling strain is independent of the longitudinal dimension and 
depends only weakly on the width of the GNR. 

It is important to note that higher e x (especially for the smaller system with small n x ) results in instabilities 
which makes the GNRs partially crumple (see Figure 5). This is due to the sp 2 bond breaking at the nonuniform 
deformed (crumpled) region of edges. Red horizontal arrow in Fig. 4(b) indicates the instability point for n x = 3. 
To find the instabilities and the buckling thresholds, we calculated the variation of the boundary forces versus 
e x . The force on the left hand side (LHS) edge is shown in Fig. 4(c). The forces on the LHS edge decrease 
to zero before the buckling thresholds which mean that GNRs at those points have an optimum length where 
the system does not feel any external forces (Fig. 5(a)). We will return to this point later. In Fig. 4(c) the 
vertical dashed line indicates the buckling threshold and the arrow gives the first sudden changes in the force 
(for three systems with n x = 2,4,6) which indicates an instability in the shape of the GNR. 

As we mentioned above, Fig. 5 shows three snap shots of a GNR (having n x = 2 size) at three different 
strains. For Fig. 5(a) strain is e = e m (strain at optimum length). Two other snap shots are GNR before 
(Fig. 5(b)) and after (Fig. 5(c)) the first instability point in Figs. 4(b,c). The strain (and compression steps) 
in Fig. 5(c) is larger than in Fig. 5(b). Both strains in Figs. 5(b,c) are larger than the critical strain (e c ). 
Brenner's potential is not responsible for the occurred non-hexagonal structures (bond breaking effects) in 
some nonuniform deformed regions of Fig. 5(c). The reason is that in common covalent potentials, people use 
a drastic reduction in cut-off distance (here ~2A) which has resulted in a good description of the material 
before fracture or bond breaking. But, it leads to an overestimation of critical loads and shear stresses in 
fracture mechanics and tribology, where bond breaking occur [15]. Here we do not study the bond breaking 
situation or fracture mechanism that can occur in GNRs by a continues compression beyond the buckling 
state. Modification to Brenner's potential have been proposed in order to include fracture mechanisms and 
bond breaking situations [16]. The idea to those modifications is to find nearest neighbors by a criterion other 
than distance [16] and this by employing empirical screening functions as introduced in ref [17]. 

3.2.1 Static buckling deformations and the stability of buckled GNRs 

In the last part of this paper, we calculate the change of the free energy of laterally supported GNRs subjected 
to axial strain. Due to the application of an external force on the boundary, an equilibrium approach is no 
longer applicable and a non-equilibrium MD is needed. Independent of the path and for a finite evolution rate, 
Jarzynski found an equality between the difference of the free energy and the total work done on the system 
(W) during a non-equilibrium evolution [7] 

AF = -/T 1 ln(exp(-/?W0), (7) 
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Figure 5: (Color online) GNR with n x — 2 size at its optimum length (a), and before (b) and after 
(c) the instability. 
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where f3 — 1/ksT. The averaging is done over the realization of the switching process between the initial and 
the final state. Equation (7) makes a connection between the difference of the equilibrium free energy and the 
non-equilibrium work. 

Using Eq. (7) we calculated the changes of the free energy when compressing the GNR of size n x = 10 
and plot the results in Fig. 6(a). The inset shows the total work done on the system for 10 simulations with 
different initial conditions. Comparing the curve for (W) and AF shows that the difference of the free energies 
are smaller than the total work which agree with (W) > AF as can be derived from Eq. (7) [7]. Fig. 6(b) shows 
the change in W for systems having different values for n x . The minimum in the free energy curve corresponds 
to an equilibrium length (also to an amount of strain (e m ) where there is zero force on the boundaries at 
x = 0,x = a (see Fig. 4(b)). 

Notice that our non-compressed GNRs (in the beginning of the simulations) are flat honeycomb lattice 
structures which are not in a thermo- mechanically equilibrium state at finite temperature. Therefore the free 
energy of this state should be higher than the equilibrium state. It is well known that at finite temperature the 
equilibrium state of suspended graphene is not exactly a flat sheet and some intrinsic ripples are present [12]. 
At the optimum length our suspended GNRs are rippled (Fig. 5(a)) and the system is in the equilibrium 
state. Figure 6(a) shows the variation of the total free energy versus applied strain. As we see from Fig. 6(a) 
the rippled state (minimum point in free energy curve) has a lower free energy with respect to the initial 
non-compressed GNRs, (AF =-2.27 eV). The inset of Fig. 6(a) depicts the corresponding change in the total 
performed work on the system for 10 simulations with different initial conditions. On the other hand, since 
the optimum length of the suspended GNR is less than the initial non-compressed length (a), we may conclude 
that at finite temperature the projected 2D-area (a(l — e m ) x b) of the GNR is less than the area of a flat GNR 
(a x b). Surprisingly, as we see from Fig. 6(b) the boundary strain at the optimum length e m = 0.48% is size 
independent. The vertical dashed lines in Figs. 6(a,b) indicate the buckling strain, i.e. e c and the strain at the 
optimum length, i.e. e m respectively. Free energy in the buckling point is less than the free energy of the initial 
non-compressed system but it is higher than the free energy for the rippled state (minimum value). Difference 
between the free energy of the rippled state and the buckled state is -1.87 eV. For the n x = 10 system we 
stopped the compression after the buckling threshold and equilibrated the system during a very long time and 
found static and stable sine wave deformations in the GNR. 

4 Conclusions 

Deformations in the graphene nano-ribbons subjected to axial boundary compression are static sine waves 
with different number of nodal lines, depending on the length of the GNRs. The deformations predicted from 
elasticity theory for the buckled rod and rectangular plate are similar to those obtained for the buckled GNRs 
in the case of free boundary condition and also for the laterally supported boundary condition. However, the 
critical force and flexural rigidity of GNRs are larger than predicted from elasticity theory. We found a linear 
relation for the stress-strain curve for small strains (i.e. before the buckling threshold). The buckling strain 
(0.86%) and the strain caused by the equilibrium length (0.48%) are independent of the longitudinal size of 
the system and they depend weakly on the width of GNRs. From the free energy of the GNRs at the buckling 
threshold, we found that they are thermodynamically more stable than those before compression, i.e. flat GNRs. 
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the system for 10 simulations with different initial conditions, (b) Total work done on the system for 
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